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The spatio-temporal dynamics of an excitable medium with 
multiple spiral defects is shown to vary smoothly with system 
size from short-lived transients for small systems to exten- 
sive chaos for large systems. A comparison of the Lyapunov 
dimension density with the average spiral defect density sug- 
gests an average dimension per spiral defect varying between 
three and seven. We discuss some implications of these results 
for experimental studies of excitable media. 
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Much research in the nonequilibrium physics of ex- 
citable media has been motivated by the observation of 
dynamical states containing defects, i. e. spiral waves 
in two space dimensions or spiral filaments in three di- 
mensions 0. Experimental studies in surface oxidation 
experiments || and in fibrillating hearts || suggest that 
many such defects may coexist in dynamically complex 
states. Although similar states have been reproduced in 
computer simulations , there has not yet been careful 
quantitative analysis of whether the long-time dynamics 
of such media can be chaotic and, if so, how the properties 
of this chaos may be related to the statistics of defects, 
to the size of the medium, and to intrinsic medium pa- 
rameters. Detailed analysis of mathematical models of 
excitable media may thus provide new insights in how 
to analyze spatially-extended excitable media, possibly 
including fibrillating cardiac tissue. 

In this Letter, we numerically study a two-dimensional 
model of a homogeneous excitable medium with an em- 
phasis on determining when spatiotemporal chaos occurs, 
and on quantitatively analyzing basic time and length 
scales of observed chaotic states. We use a model in- 
troduced by Bar et al. as a reduced description of 
carbon monoxide oxidation on a surface [Q, because of 
its numerical simplicity and because of prior work sug- 
gesting the existence of chaos B. Extending some recent 
work by other researchers (J,|5|,|9[ , we find that the dy- 
namics is strongly dependent on the system size L. For 
small L, all initial conditions studied rapidly decay to 
an asymptotic constant or periodic state. As the system 
size increases, however, the volume of the set of initial 
conditions leading to sustained, non-periodic dynamics 
increases smoothly, and we discuss the transition from 
periodic to non-periodic dynamics with increasing system 
size. The non-periodic dynamics sustained in sufficiently 



large systems are statistically stationary, and we compute 
Lyapunov exponents and dimensions, defect statistics, 
and two-point correlation lengths to characterize these 
states. These different statistics are compared, both to 
test previous conjectures about the relationship of these 
correlation lengths Jl(| and to evaluate the complexity of 
the defects. Our results indicate that a two-dimensional 
excitable medium of moderate size with few defects on 
average can already sustain extensive, high-dimensional 
chaotic dynamics, a fact with important implications for 
control of excitable media by small parameter perturba- 
tions [pd) . In the following, we explain the model, sum- 
marize our calculations, and discuss the implications of 
our results. 

The Bar model describes the interaction of an activa- 
tor field u(t, x, y) with an inhibitor field v{t, x, y) via the 
partial differential equations 
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which we solve numerically in a square domain of side L 
with either biperiodic (BP) or no-flux (NF) boundary 
conditions on the field u(t,x,y). The function f(u) has 
the form 
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if u< 1/3, 
if 1/3 < u < 1, 
if u > 1, 
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so that the production of the inhibitor v is "delayed" 
until u exceeds 1/3. The nonlinear form Eq. || leads to 
three fixed points, one stable and two unstable; the larger 
unstable fixed point (u*,v*) (which does not appear in 
the widely-used Fitzhugh-Nagumo model) seems neces- 
sary for the occurrence of spatiotemporal chaos. The 
parameter e in Eq. ^ determines the ratio of time scales 
of the fast field u and slow field v and is the key bifurca- 
tion parameter in this paper. The positive parameters a 
and b were fixed at the values a = 0.84 and b = 0.07 to 
take advantage of substantial earlier work using these val- 
ues Spiral solutions are then known empirically to 
be unstable when e exceeds a critical value e c ~ 0.069 [0. 
The mechanism of this instability, meander of the spiral 
core into a branch of the spiral, is apparently unique 
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to models with "delayed" inhibitor production like that 
given by Eq. (||). In particular, this is not the mecha- 
nism of breakup observed in models of cardiac tissue |&H . 
Breakup leads to long-lived, complicated dynamics for 
certain initial conditions when e > e c ; a snapshot of such 
a disordered nonperiodic state with 31 spiral defects is 
shown in Figure 1. 

Our calculations involved integrating Eq. (1), calculat- 
ing the Lyapunov spectrum of the numerical trajectory, 
and counting the number of spiral defects at successive 
times. For both kinds of boundary conditions, Eq. (1) 
was solved numerically by first introducing second-order 
accurate finite difference approximations for the spatial 
derivatives on a uniform square mesh of spacing A a; and 
then using an algorithm proposed by Barkley [fl2[ . For 
the calculations reported below, we used a spatial grid 
size Ax = 0.50 and time step At = 0.05e. The spec- 
trum of Lyapunov exponents Xi and the Lyapunov fractal 
dimension D were calculated by well-known algorithms 
based on linear variational equations |l5| l that were inte- 
grated by a forward-Euler algorithm with the same grid 
and time step. The time step chosen, At = 0.05e, was 
much smaller than that required by integration of only 
Eq. (1), but was found to be necessary to compute Lya- 
punov exponents accurately to within a few percent p[ . 
For given boundary conditions and initial data, Eq. (1) 
was integrated for 2000 time units to allow a statisti- 
cally stationary state to be obtained, and then the full 
system with variational equations was integrated for an 
additional 1000 time units (~ 200 spiral periods), during 
which statistics were calculated. 

To study the dependence of the dynamics on initial 
conditions, we integrated Eq. (1) from each of 100 ini- 
tial conditions generated by distributing the field val- 
ues uniformly (at each grid point) in the ranges u 6 
[0.8u*,1.2u*],v S [0.8u*,1.2«*]. This procedure was re- 
peated for both boundary conditions and for square sys- 
tems of side length L varying from 5 to 40. For all initial 
conditions, the dynamics was short-lived in small sys- 
tems (L < 15 (for NF boundary condition) or L < 8 
(for BP) ), decaying in less than 100 time units to either 
the stable uniform state or to a plane-wave state (only 
in the case of biperiodic boundary conditions). Suffi- 
ciently large systems {L > 35) sustained dynamics for 
at least 3000 time units. The fraction / of initial con- 
ditions which led to non-periodic dynamics sustained for 
a time T (either 100 or 1000 time units) is shown in 
Figure || as a function of system size for both boundary 
conditions. For biperiodic boundary conditions (Figure 
2(a)), the curve is independent of the cutoff time T, indi- 
cating that transients either die quickly or are sustained 
indefinitely (more than 50,000 time units). For no- flux 
boundary conditions, all initial conditions studied even- 
tually decayed to a stationary state; the mean and me- 
dian transient times both scale exponentially with sys- 
tem size pS], much like the supertransient behavior ob- 



served previously in one-dimensional systems jl5[ . These 
results suggest that excitable media of intermediate size 
may have an appreciable basin of attraction both for non- 
periodic dynamics and for periodic or constant dynam- 
ics. Whether this accounts for the observation that fib- 
rillation sometimes occurs in hearts of intermediate size 
remains unclear, both because of the differing breakup 
mechanisms and because of the effect of the third dimen- 
sion in heart tissue [p|,^6[. 

The fact that a given state was transient was revealed 
only by an eventual abrupt change to the uniform state; 
the dynamics of the transient itself was found to be sta- 
tistically stationary. For parameter values e > e c and 
for system sizes L > 25, these statistically stationary 
states were found to be high-dimensional (D > 20) and 
extensively chaotic [jlTj as shown in Figure 3 by a lin- 
ear dependence of D on L 2 . From the asymptotic slopes 
of the curve in Figure 3, an intensive dimension den- 
sity 5 — lim^2^ oc dD/d(L 2 ) was obtained and then re- 
expressed as a dimension correlation length £g = S^ 1 ^ 
for a d = 2 dimensional domain (lj] . To test a spec- 
ulation of Bayly et. al. |[o| that knowledge of the 
experimentally-accessible two-point correlation length £2 
might provide knowledge of the dynamical length £s for a 
chaotic state of spiral defects, we computed £2 and £5 for 
several values of the parameter e. For each e value stud- 
ied, the two-point correlation function C(r) had a sim- 
ilar monotonically-decreasing but non-exponential form 
so we estimated £2 by the position of the first zero cross- 
ing of C(r) iQ. As shown in Figure 4(a), the two lengths 
agree within a factor of 1.5 or better but have opposing 
trends as e increases (from 0.07 to 0.12), with £5 decreas- 
ing and £2 increasing slightly. It is unclear from this data 
whether an estimate of £5 can, in general, be obtained by 
measuring £2- In any case, further analysis of more phys- 
iologically accurate models will be needed to relate £g 
and £2 for the heart data of Bayly et. al. [|Io| . 

We also explored whether the fractal dimension D of 
the chaotic states was related to the statistics of the num- 
ber N(t) of spiral defects, e.g., to its time average (TV). 
The spirals were counted at successive times by locat- 
ing their cores, which occur at those points (x,y) in 
the medium where the fields (u, v) take on the values 
(u*,v*) 0. It has been shown previously that N(t) is 
constant for e < e c , in which case the time average (N) is 
fixed by the choice of initial condition Q . We found that 
the mean (TV) scaled extensively with system area L 2 for 
both boundary conditions considered, so that the average 
defect density n = (N) / L 2 was independent of system 
size ||]. The ratio d = D/(N) of the Lyapunov dimen- 
sion to the mean number of defects therefore defines an 
intensive quantity that measures the number of dynam- 
ical degrees of freedom associated with each defect on 
average. Because the extensive quantities D and (iV) are 
of the form aL 2 + /3 rather than simply aL 2 (where a 
and (5 are constants), the ratio d asymptotes slowly to a 
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constant value close to S/n. We studied the dependence 
of D/(N) on area L 2 for two different values of e, and 
found that we could estimate the ratio d accurately us- 
ing a single system of side length L = 40 with periodic 
boundary conditions. For this system size, d increases 
smoothly with e from less than 4 for e ss e c to nearly 7 
for e > 0.095 at which point it becomes approximately 
constant (see Figure 4(b)). Thus a fixed number of de- 
grees of freedom can not, in this manner, be associated 
with each spiral defect in an excitable medium. 

In summary, for a particular model of a two- 
dimensional excitable medium Q^], we have demon- 
strated by numerical calculations a smooth transition 
from short-lived transient dynamics to extensive, high- 
dimensional chaotic dynamics with increasing system 
size L for both no-flux and biperiodic boundary condi- 
tions. Small systems never exhibited sustained chaotic 
dynamics, but non-periodic dynamics in sufficiently large 
systems were found to be statistically stationary on such 
long time scales and for such a large fraction of random 
initial conditions that Lyapunov spectra of the dynamics 
converged well to a value independent of the initial con- 
dition. Although previous results based on time series 
allowed computation of one Lyapunov exponent Jl8| , p^| , 
we computed enough Lyapunov exponents to determine 
the Lyapunov dimension D in an excitable medium with 
many spiral defects. We found no clear relation between 
the dimension correlation length £5 and the widely used 
two-point length £2; however, the numerical similarity of 
these lengths supports a previous conjecture that £g w £2 
in data taken on fibrillating hearts |l0| . The mean Lya- 
punov dimension per defect of 3 to 7 suggests that defects 
are more complicated than in the defect-turbulent regime 
of the complex Ginzburg-Landau equation |20| ] , and that 
the dynamics of excitable media with even a few defects 
may be quite high-dimensional. This high-dimensionality 
in turn suggests that it will be difficult to stabilize such 
states by small variations of parameters |Tl[] , and may 
explain why some previous attempts to analyze the dy- 
namics of fibrillation with low-dimensional time series 
embedding techniques have been inconclusive [2]}| . The 
unusual nature of the spiral- wave breakup in this medium 
leaves to future work the important question of whether 
the results obtained here apply to other excitable media, 
including fibrillating ventricles. 
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FIG. 1. Density plot at time t = 500 of the slow field 
v(t,x,y) for a spatiotemporal chaotic state with 31 spiral de- 
fects present. Dark and light regions correspond to values less 
or greater than the value v* — 0.484 corresponding to the un- 
stable fixed point; the field values span the range v € [0, a — b\. 
Parameter values were e = 0.074, a = 0.84, b = 0.07, L = 50, 
Ax = 0.5 and At = 0.0037. 



FIG. 2. Fraction / of 100 random initial conditions still 
exhibiting non-periodic dynamics after a given time, (a) For 
biperiodic boundary conditions with cutoff time T np = 100 
(circles), 1000 (squares), or any larger value, the transition 
has the same form, with systems of side length L > 25 
nearly always exhibiting sustained dynamics, (b) For no-flux 
boundary conditions with cutoff time T np — 100 (circles) and 
Tnp = 1000 (squares), we see that the median transient time 
depends on the cutoff. Comparison of these graphs shows 
that dynamics are substantially less likely to be sustained for 
a given time with no-flux boundary conditions. The parame- 
ters used were the same as in Figure 1. 



FIG. 3. Lyapunov dimension D versus system area A — L 
of Eq. (1) for the parameter values of Figure [[]. Extensive 
(linear) scaling is found for two different boundary conditions, 
no-flux (squares) and periodic (circles). Data for L < 25 
did not exist since all initial conditions decayed quickly to 
the uniform state. The dimension extrapolates to zero for a 
positive system size, so the ratio D/L 2 of the dimension to 
system area asymptotes slowly to the dimension density S. 



FIG. 4. (a) Dimension correlation length £a (circles) and 
two-point correlation length £2 (squares) for different e values, 
(b) Degrees of freedom per mean defect, D/(N) as a function 
of e. This ratio increases steadily with e above the transition 
to chaos at e — e c (marked by the arrow), and varies little in 
the region e > 0.095. 
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